# generate a realization of heterogenous LGCP model 
# 
# Author: Guochun Shen
# Data:   2011-12-18
# Project:spatial statistic
# Email:  shenguochun@gmail.com
###############################################################################



rHeteLGCP=function(en.filter,sigma2,alpha,N,plotdim,nu=0.5){
	sigma2true=sigma2
	covrdim=dim(en.filter$v)
	xcell=plotdim[1]/covrdim[2]
	ycell=plotdim[2]/covrdim[1]
	mu=log(N/sum(exp(en.filter$v+0.5*sigma2true)*xcell*ycell))
	xcol=seq(xcell/2,plotdim[1],xcell)
	yrow=seq(ycell/2,plotdim[2],ycell)
	retry=TRUE
	r=seq(0,100,1)
	while(retry){
		if(nu!=Inf){
			Y <- GaussRF(x=xcol, y=yrow, model="matern", grid=TRUE,
					param=c(mean=0.0, variance=sigma2true, nugget=0.0, scale=alpha,nu=nu))
		}else{
			Y <- GaussRF(x=xcol, y=yrow, model="gauss", grid=TRUE,
					param=c(mean=0.0, variance=sigma2true, nugget=0.0, scale=alpha/sqrt(0.5)))
		}
		Yim <- as.im(list(x = xcol, y = yrow, z = Y))
		#mu2=log(N/sum(exp(Yim$v)*xcell*ycell))
		#error=checksimpp(Yim,r,c(sigma2,alpha),nu,mu2)
		#if(error<2){
		retry=FALSE
		#}
	}
	
	Lambda=en.filter
	Lambda$v=exp(mu+en.filter$v+Yim$v)
	#simulate inhomogeneous Poisson process with intensity function given by Lambda.
	X=rpoispp(Lambda)
	return(X)
}

